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ABSTRACT 

We use the Delaunay Tessellation Field Estimator (DTFE) to study the one-point density 
distribution functions of the Millennium (MS) and Millennium-II (MS-II) simulations. The 
DTFE technique is based directly on the particle positions, without requiring any type of 
smoothing or analysis grid, thereby providing high sensitivity to all non-linear structures re- 
solved by the simulations. In order to identify the detailed origin of the shape of the one-point 
density probability distribution function (PDF), we decompose the simulation particles ac- 
cording to the mass of their host FoF halos, and examine the contributions of different halo 
mass ranges to the global density PDF. We model the one-point distribution of the FoF halos 
in each halo mass bin with a set of Monte Carlo realizations of idealized NEW dark matter 
halos, finding that this reproduces the measurements from the N-body simulations reasonably 
well, except for a small excess present in simulation results. This excess increases with in- 
creasing halo mass. We show that its origin lies in substructure, which becomes progressively 
more abundant and better resolved in more massive dark matter halos. We demonstrate that 
the high density tail of the one-point distribution function in less massive halos is severely 
affected by the gravitational softening length and the mass resolution. In particular, we find 
these two parameters to be more important for an accurate measurement of the density PDF 
than the simulated volume. Combining our results from individual halo mass bins we find that 
the part of the one-point density PDF originating from collapsed halos can nevertheless be 
quite well described by a simple superposition of a set of NEW halos with the expected cos- 
mological abundance over the resolved mass range. The transition region to the low-density 
unbound material is however not well captured by such an analytic halo model. 
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1 INTRODUCTION 

One of the most important questions in cosmology is to understand 
the formation of large-scale structures in the Universe. In the stan- 
dard ACDM model, the energy density of today's Universe is dom- 
inated by non-baryonic cold dark matter (~ 23%) and dark energy 
73%), whereas only ~ 4% is comprised of all the mass and 
energy associated with planets, stars, galaxies, clusters, gas, dust 
and electromagnetic radiation. The study of large scale structures 
is hence primarily the study of the distribution of galaxies and the 
underlying dark matter. 

The distribution of galaxies can be well studied w ith large 
gala xy redshift surveys (e.g. SDSS lStoughton et aLlJiool , or 2dF- 
GRS IColless et al.l bOOlh ) which provide maps of the distribution 
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of galaxies over large volumes in the (mostly) nearby Universe. On 
the other hand, there is no direct detection of dark matter particles 
yet and its existence is mostly inferred from indirect observations 
such as gravitational tensing, making its observational study much 
harder. In the current paradigm of structure formation, structures 
form via gravitational instability amplifying tiny density fluctua- 
tions generated by some process in the early Universe. These initial 
density fluctuations are often assumed to form a Gaussian random 
field. Dark matter first aggregates hierarchically into dark matter 
halos and galaxies form later in their centres by the cooling and 
condensation of baryons (White & Rees 1978). As the dark matter 
halos form near peaks of the initial density field, the distribution of 
dark matter and galaxies on large scales is larg ely determined by the 
statistics of these peaks jSardeen et alj|l98q) . On small scales, the 
physics of galaxy formation complicates this picture considerably. 
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leading to non-linear and stochastic biasing between the distribu- 
tions of dark matter and galaxies. 

Studies of the cosmic density field expected in cold dark mat- 
ter cosmologies are often based on s imple and approximate analyt- 
ical models such as the halo model dCoorav & Sheflill2002 l). How- 
ever, detailed studies of the non-linear cosmic density field need to 
rely on N-body simulations, which do not need to make simplify- 
ing assumptions about the abundance and structure of halos. In the 
present work, we study the non-linear density fields predicted by 
high-resolution dark matter simulations, particularly the one-point 
probability distribution function of the dark matter density field. We 
measure this function far into the nonlinear regime and compare the 
results to the halo model. 

The output of N-body simulations provides the phase space 
distribution of dark matter particles. Reconstructing the under- 
lying continuous density field represented by the discrete set of 
"macro-particles" used by the numerical scheme requires one to 
define an appropriate density reconstruction scheme. For consis- 
tency, we demand that the total mass contained in the recon- 
structed continuous density field has to be exactly equal to the 
total mass represented by the discrete set of particles. There ex- 
ist various techniques in the literature for density reconstruc- 
tion from a given set of points which fulfill t hese requirements 
(Hocknev & Eastw ooj Il98ll : Isilvermaj 1 19861 : iMonaghan 1992; 
IXscasibar & Binnevll2005b 7The most widely used approach is to 
convolve the point data with some filtering function (or simply 
'kernel'), yielding a continuous map. Conventionally, the filtering 
function has a fixed shape and fixed size (for example when bin- 
ning particles on a regular grid by CIC or TSC mass assignment), 
but this fixed smoothing technique has the serious disadvantage that 
the smoothing length is not adjusted to clustering of the particle dis- 
tribution. So when a small smoothing length is employed in order to 
achieve great resolving power in high density regions like filaments 
and clusters, these structure are well recovered but underdense re- 
gions like voids are severely affected by shot noise. Conversely, 
if one wishes to obtain a reasonable reconstruction of low-density 
voids by using a larger smoothing length, the filaments and clusters 
are oversmoothed, limiting the amount of information that can be 
extracted from those regions. 

A better smoothing technique is obtained by applying the SPH 
(smoothed particle hydrodynamics) approach, where one employs 
an adaptive kernel which adjusts itself according to the varying 
sampling density. For example, the size of the (compact) kernel 
can be set to the distance of the n-\h nearest neighbor, where the 
value of 71 is a user-specified parameter. In both of these meth- 
ods, the smoothing kernel has to be specified by the user. Common 
choices consist of spherically symmetric kernels, for example a 
simple Gaussian. However, the fact that the geometry of the kernel 
is prescribed irrespective of anisotropics present in local non-linear 
structures (e.g. filaments and sheets) may introduce spurious topo- 
logical signatures characteristic of the kernel. Ideally, one would 
like to allow the point distribution to decide for itself what kernel 
shape and size yields the most faithful reconstruction of the local 
density field. An ideal candidate for this strategy is the Voronoi 
tessellation and/or its t opological du al, the Delaunay tesse llation 
("van de Wevgaert' 1994; Okabe et al.l l200(t IPelupessv et"Zll2003l : 
[van de Weygaert 2007; Schaap 2007). The density estimators based 
on these tessellations have several advantages over the traditional 
smoothing techniques which we discuss in the next section. These 
advantages of DTFE over traditional smoothing techniques for den- 
sity reconstruction have made it an increasingly popular choice 
in recent years (e.g. lAragon-Calvo et aLll2007llzhang et al.ll201Cl 



|Platenetal.]|201lh . We will primarily use the Delaunay tessella- 
tion field estimator (DTFE) because it offers a parameter free re- 
construction of the density field, retaining a maximum amount of 
information about the density field and the topology of structures 
embedded in the dark matter distribution. 

The Millennium Simulation (hereafter MS) ( ISpringel et al.l 

I2OO5I) is still one of the largest high-resolution simulations of 
the growth of dark matter structures. It followed the evolution 
of 10 billion dark matter particles in a 500^~^Mpc comoving 
box with an individual particle mass of 8.6 1 X 10^ h-^ Mp). The 
Mille nnium-II simulation (hereafter MS-II ; lBovlan-Kolchin et al.l 
I2OO9I) simulated a 100 /i ^Mpc box using the same number 
of particles, thereby offerin g 125 times better mass resolution. 
iBovlan-Kolchin et alj f2009l) studied the formation and statistics of 
dark matter halos in the MS-II simulation. By comparing their re- 
sults with the MS simulation they found excellent convergence in 
the basic dark matter halo statistics, making these two simulations 
ideally suited for a study of the dark matter density field over an 
unprecedented range of scales. 

The one point distribution function of the cosmological den- 
sity field is one of most fundamental quantities characterizing sta- 
tistical properties of the matter distribution in the Universe. In the 
current paradigm of structure formation, the present day large scale 
structures grew from primordial density fluctuations with Gaus- 
sian statistics. The one point distribution of today's density field 
is however far from Gaussian as a result of gravitational evolu- 
tion. In the mildly non-linear regime, it is known that the one point 
distribution of the dark matter density field obtained from N-body 
simu lations is reas onably well described by a log-normal distribu- 
tion dColesfe JonesHl99ll; iKofman et aljri994l ; TKavo etal]|200ll ; 
iTaruva et al.l2003h . but this approximation eventually breaks down 
in the highly non-linear regime. 

In this paper, we study the one point distribution of the dark 
matter density fields in the MS and MS-II using DTFE and try 
to interpret the results in the simple picture provided by the halo 
model. The dark matter halos are the densest sites in the cosmic 
mass density field, and approximately ~ 50% (49.6% in the MS 
and > 60% in the MS-II at 2 = 0, with a 20 particle limit) of the 
mass is bound in resolved halos. Note that in the halo model, all of 
the mass in the Universe is assumed to be part of a dark matter halo 
of some mass. The density profiles of simu lated CDM halos ar e 
well described by the universal NEW profile dNavarro et al.lll99^ . 
with a shape approximately independent of mass, the amplitude 
of initial density fluctuations an d cosmology dNavarro et alj|l997l ; 
ICole&Lacevll 19961 : |jin3l200(]|) . The concentration varies weekly 
with halo formation time. In the halo model, one tries to represent 
the underlying dark matter distribution as a superposition of a set 
of NEW halos with abundance and clustering modelled with sim- 
plistic models or analytic fits to N-body results. 

N-body simulations compute a periodic model universe of fi- 
nite size and finite mass resolution. This also requires a soften- 
ing length below which the gravitational interaction is suppressed 
to avoid singularities in orbit integrations and unphysical particle 
scattering. These numerical limitations are expected to influence 
the ability of the simulation to resolve very high and very low 
density regions, and consequently affect the tails of the one-point 
distribution. The MS and MS-II use different simulation volumes, 
mass resolutions and softening lengths, allowing us to study the 
importance of these effects in shaping the tails of the one-point 
distributions. We note that the use of a user defined kernel for es- 
timating densities like in SPH introduces a smoothing scale and 
a corresponding resolution element which will typically have an 
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Figure 1. The probability distribution functions of VTFE and DTFE recon- 
structed density fields for a Poisson point process with lO'' particles in 3D. 
The one-point distribution function in each bin is weighted by the bin width 
and the volume weighted density in that bin. An analytic approximation for 
the VTFE reconstructed density field given by Eqn. (|3] is shown as a solid 
line. The coiTesponding best fit values a, b, and c for VTFE and DTFE are 
listed in Table 1 . 

additional effect on the tails of the distribution. As DTFE is self- 
adaptive without a free parameter, this type of effect is expected to 
be less important for this scheme than other numerical limitations 
due to finite volume, finite mass resolution and gravitational soft- 
ening, additional effect on the tails of the distribution. As DTFE is 
self-adaptive without a free parameter, this effect is less important 
than other numerical effects for finite volume, finite mass resolu- 
tion and gravitational softening. 

This paper is structured as follows. In Section |2] we outline 
the Delaunay Field Tessellation Estimator and discuss our method- 
ology. We then present results for toy halo models and the MS and 
MS-II simulations in Section [S] Finally we conclude and summa- 
rize our findings in Section|4] 



2 DENSITY ESTIMATION WITH THE DTFE 

In mathematics and computational geometry, the Delaunay tessel- 
lation for a set of points is the uniquely defined volume-covering 
tessellation with mutually disjoint tetrahedra, in which no circum- 
sphe re of any tetrahedron contains on e of the points in its inte- 
rior (lDeIaunavll934IOkabe et al.l2000l) . Connecting the centers of 
the circumscribed spheres of neighboring Delaunay tetrahedra pro- 
duces the Voronoi tessellation of the point set, which is the topolog- 
ical dual of the Delaunay tessellation. The Voronoi tessellation is a 
division of space into non overlapping convex regions where each 
region is uniquely assigned to one of the sampling points. All the 
points in these convex regions are closer to its defining sampling 
point than to any other sampling point. 

Based on the geometric constructions of these tessellations, 
different density reconstruction schemes can be constructed. For 
example, the density at each sampling point in the VTFE (Voronoi 
Tessellation Field Estimator) is simply defined as pi = rrii/Vi, 
where is the mass of the i-th sampling point and Vi is the 
volume of the corresponding Voronoi cell. This method assumes 
that the mass of each particle is uniformly distributed inside each 
Voronoi cell, keeping the density constant inside each cell. The 
product of density and volume of all the Voronoi cells trivially re- 
turns exactly the total mass of all the sampling points. But an im- 



Table 1. The number of points N, and the fitting parameters a, b and c 
assummg a functional form /(p) = a p-'' e^^^P for the PDF of VTFE 
and DTFE reconstructed density fields for a 3D Poisson point processes. 
The mean is 1 for both PDFs. The variances of the PDFs of VTFE and 
DTFE are 0.22 and 0.24 respectively. 



portant deficiancy of this density reconstruction is that the density 
field is discontinuous at the Voronoi cell boundaries. 

An improved density estimator that addresses this deficit is 
the Delaunay Tessellation Field Estimator (DTFE), which is based 
on a Delaunay tessellation of the sampling points, as proposed by 
ISchaap & van de Weygaertl feOOd) . Here the density at each sam- 
pling point is defined as pi = 4 m; /Wi where Wi is the volume 
of the contiguous Delaunay region around the point (composed of 
all the tetrahedra that have the point as one of their vertices). The 
sum Wi — is the sum of the volumes of all Delaunay 

tetrahedra that share point i as one of their vertices. The multipli- 
cation by 4 accounts for the fact that each Delaunay tetrahedron is 
contributing to the contiguous Delaunay region of four points. The 
DTFE density estimation scheme assumes that the density field in- 
side each tetrahedron varies in a linear fashion. The gradient of the 
density is assumed to be constant within each tetrahedron and can 
then be computed using the density values at the four vertices of the 
tetrahedron. One can then easily find the density at any other loca- 
tion inside the tetrahedron using tri-linear interpolation. This cre- 
ates a continuous, volume covering, piece-wise linear density field. 
It is easy to verify that the volume integral of the DTFE density 
field reproduces the sum of the particle masses exactly. The most 
important advantage of DTFE over conventional methods is that 
the density estimates in this method do not rely on any additional 
parameter. The DTFE kernel not only adapts to the local density as 
in the case of SPH but also to the local geometry of the distribu- 
tion. We emplo y the tessellation engine of the parallel AREPO code 
jSpringelllOloh to construct the Delaunay mesh. 

To construct and store the Delaunay mesh, AREPO uses in- 
dices to refer to 4 vertices and 4 adjacent tetrahedra of each tetra- 
hedron which requires at least 32 bytes of memory per tetrahedron 
on 64-bit machines (plus 4 bytes for an auxiliary variable in prac- 
tice), provided the number of points per distributed memory region 
is kept low enough to allow the use of 32-bit integers (which is 
the case in practice). For a r andom point set there a re on average 
~ 6.77 tetrahedra per point (Ivan deWevgaert|[T994l) , implying at 
least ~ 244 bytes of memory per point for storing the mesh tetra- 
hedra. Another 20 bytes per point are required to hold the parti- 
cle coordinates (if stored in single precision) and a unique parti- 
cle ID. In practice, additional memory is needed for a search tree 
(in order to validate individual Delaunay tetrahedra by efficiently 
searching for points inside the circumsphere) and for 'ghost' cells 
that mesh the different tessellation patches together across proces- 
sor domains. There are more than 10 billion particles in both MS 
and in MS-II. This requires us to build a Delaunay tessellation com- 
posed of more than 70 billion tetrahedra for each of these simula- 
tions. We used the ODIN machine at the Computing Center of the 
Max Planck Society, Garching, to perform the mesh constructions 
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Figure 2. Distributions of p/pNFW for all the particles within r2oo in an 
NFW halo with parameters A/200 = 10^ M0, Af200 = 10®, r200 = 
50 kpc and c = 5 (black symbols). The distributions for particles residing 
within certain radial ranges are separately shown, as indicated in the fig- 
ure. Note that the distributions shown in this plot ai'e simply the histograms, 
and their values in each bin ai'e not weighted by the bin width or the volume 
weighted density as in Figure[T] The distributions in different radial bins are 
all of the same shape and are consistent with what we obtained for the Pois- 
son sampling of a uniform distribution, even though the sampling densities 
vary strongly across the different radial bins. 
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for MS and MS-II, using 512 cores and ~ 7.5 Terabytes of memory 
in total. In both cases, the mesh construction took about 20 minutes 
of wall-clock time. 



3 RESULTS 

In cosmology, we quite often encounter Poisson sampling of an un- 
derlying density field. For example, the galaxies in a redshift survey 
or the particles in an N-body simulation can be considered as Pois- 
son samples for certain density ranges. Therefore it is important to 
understand the impact of the Poisson sampling noise on the statis- 
tics of the reconstructed density distribution. Both the VTFE and 
the DTFE reconstruct the density field with an adaptive spatial res- 
olution from a discrete set of data points. The high sensitivity of 
these density estimators to the variation of local density and geom- 
etry makes them presumably particularly sensitive to the presence 
of shot noise. We therefore first examine the statistical properties 
of Poisson sampling noise as seen by these density estimators. Fur- 
ther, the datasets in cosmology also often involve Poisson sampling 
of highly inhomogeneous distributions for example from N-body 
simulations and galaxy surveys. It is also important to test whether 
the one point distribution for uniform Poisson sampling of a ho- 
mogeneous distribution also describes the noise caused by Poisson 
sampling an inhomogeneous distribution. We address this question 
with Monte Carlo simulations of NFW halos. 

We start this section by characterising the performance of the 
DTFE for a Poisson sample (in Subsection 13. Il l and for Monte 
Carlo realisations of NFW halos (in Subsection 13. 2t . These test 
cases, for which the underlying continuous density field is known, 
will help us to understand and interpret the main results of the pa- 
per, the one point density distributions of the cosmic density field 
given by the DTFE applied on the Millennium simulations (in Sub- 
section [33t . We finalise this section by examining the impact in our 
results of halo substructure and ellipticity, as well as of numerical 
setup of the simulations. 
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Figure 3. The top panel compares the distributions of p/pnewC^'s) for all 
particles in a NFW halo with parameters M200 = 10® Mq, N200 = 10®, 
r2oo = 50 kpc, c = 5, and in all the 10"* halos with Af2oo = 10® M©, 
^200 = 10^, r200 = 50 kpc, c = 5. The bottom panel shows the cu- 
mulative mass as a function of density in the halo. The results show that 
the DTFE can quite reliably represent the mass profile of the entire halo 
even with a particle number as low as N200 = 10^. The overprediction 
of the cumulative mass at r < in the Af200 = 10® halo is caused by 
confiibutions from other radii due to large scatter in the DTFE densities. 

3.1 Density PDF of a 3D Poisson point process 

We begin by studying the one-point distribution function of the 
density for a Poisson process analyzed with Voronoi tessellations. 
A Poisson Voronoi tessellation results when the generating points 
of the Voronoi cells are a Poisson point sampling of a uniform field. 
In the case of one-dimensional Voronoi tessellations, one can rig- 
orously derive the probability distribution of the lengths of the seg- 
ments, which is given by 

5r(i) = 4i exp(-2i), (1) 

where x — l/{l). Here I is the length of the Voronoi cell and (l) its 
average. No analytical results are known for the size distributions 
of Poissson Voronoi cells in 2D and 3D. Empirical studies using 
Monte Carlo realizations fit the distribution of surface area or vol- 
ume of the Voronoi cells in 2D or 3D with a gamma type probability 
distribution function (Kiang 1966) approximated by 

9{^;a,) = -——x e , (2) 
r(a) 

where x = v/{v) is the size of the Voronoi cell in units of the 
average cell size, > i > 00, and a is a constant whose value de- 
pends on the dimensionality of the space. Monte Carlo experiments 
suggest o = 2, 4 and 6, for ID, 2D and 3D, respectively. 
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The probability of a random point to lie inside a Voronoi cell 
of size X is the product of g(x) and x, which in other words is 
the probability of a random point to have a density p = 1/i in 
the VTFE reconstructed density field. Following this definition, the 
one point distribution function of the VTFE reconstructed density 
field in 3D is 

Ag{p) = dg{x) xx = 388.8 p"^ e~^''^ dp, (3) 

if the hypothesis that g{x\ a) follows equation ^ with a = 6 is 
indeed correct. To verify these results, we generate a Poisson point 
process with 10® points in a cubic box 100 Mpc on a side. We 
construct the Voronoi tessellation and Delaunay tessellation of the 
points and then estimate the density field values at the sampling 
points using both VTFE and DTFE. We compute the one point dis- 
tribution of both the VTFE and DTFE reconstructed density fields. 
We use only the density estimates at the sampling points for com- 
puting the one-point distribution of the VTFE reconstructed density 
field, and determine the best-fit parameters assuming a functional 
form /(p) — ap'^eT'^^^ to describe the one-point PDF of the 
VTFE reconstructed density field. In the VTFE density reconstruc- 
tion scheme, the density at the location of the sampling points is 
defined as the inverse of the volume of its coiTesponding Voronoi 
cell weighted by its mass, whereas in the DTFE scheme the vol- 
umes of the contiguous Delaunay cells are used instead. 

In the DTFE density reconstruction, the density field inside 
each tetrahedron varies linearly, whereas in the VTFE density re- 
construction the density inside each Voronoi cell is constant. To 
account for this difference and arrive at a more appropriate com- 
parison between DTFE and VTFE, we randomly select one point 
inside each Delaunay tetrahedron through uniform sampling and 
determine its density estimate by linearly interpolating from the 
four vertices of the corresponding tetrahedron. We weight each 
such point by the volume of the corresponding tetrahedron. We also 
weight each point in VTFE by the volume of their Voronoi cells. 
As there are on a verage 6.77 Delauna y tetrahedra per point for a 
random data set jvan de WeYgaert|[l994l) . we get 6772467 density 
estimates from the DTFE reconstructed density field by randomly 
choosing one point from each tetrahedron. The one-point distribu- 
tions of the VTFE and DTFE reconstructed density fields are both 
fitted to the same functional form of equation (|3} and the best-fit 
parameters are listed in Table 1. We find that the distributions ob- 
tained for the VTFE density field are well described by equation 
(|3), and the values of our best-fit parameters approach the literature 
values as the number of particles is increased. The distribution ob- 
tained for the DTFE density field is pretty similar overall, but it is 
clearly not exactly the same. Both for the VTFE and DTFE recon- 
structions, equation (|3} noticeably undeipredicts the high density 
tail of the one-point distribution function. It is also interesting that 
the two schemes have very similar variances despite the seemingly 
larger smoothing involved in the DTFE scheme. 



3.2 Monte Carlo realizations of NFW halos 

Dark Matter halos in N-body simulation are highly inhomogeneous 
system s which we idealize as sy stems with spherical NFW density 
profile (Navarr o et aLlll996l[l997l) . 

The NFW density profile is given by 

= TT . (4) 



corresponding to a cumulative mass within radius r of 

M(t) = 1^ 47Tr'^p{r') dr' = 47Tp,5,rt [ln(l + ^) - j^^] 

(5) 

Here pc = 3Hq / (SnG) is the critical density, — r2oo/c is the 
scale radius, r2oo is the virial radius, c is the concentration of the 
halo, and 5c is the characteristic density. The virial mass of the halo 
is -A/200 = 200 Pc (47r/3) r|oo- The characteristic density 

_ 200 c3 

3 ln(l + c) -c/(l + c) ^ ' 

is related to the concentration c by the requirement that the mean 
density within r2oo should be 200 times the critical density. 

To generate our mock NFW halos, we first specify M200 and 
the particle mass. This determines the number N200 of particles 
which reside within r2oo. Next, we specify the concentration of 
the halo using the mass-co ncentration relation determined from the 
Millennium Simulation by Neto et al. (2007). We then populate the 
halo with particles using a Monte Carlo sampling technique, i.e. the 
probability to place a particle at a certain radius is made propor- 
tional to dM from equation (|5}. This is augmented with isotrop- 
ically selected angular co-ordinates 9 and (j). To avoid boundary 
effects, we extend the NFW halo radially out to 3 r2oo- 

3.2.1 One-point distribution in dijferent parts of the same halo 

We generate Monte Carlo realizations of an NFW halo with the fol- 
lowing parameters: Af2oo = 1O®M0, N200 — 10'', r2oo ~ 50kpc 
and c = 5. The halo extends out to 3 x r2oo- We put the halo at the 
center of a cubic box with a side of length 6 x r2oo and construct the 
Delaunay tessellation using all the particles in the halo. Note that 
this leaves empty regions at the comers of the box. The Delaunay 
cells near the halo boundary will then be very extended due to the 
presence of these empty regions resulting in spurious density esti- 
mates. In order to avoid these boundary effects we limit our analysis 
to the particles residing within the virial radius r2oo of the halo. We 
choose three radial bins, r/r^ = 0.01 — 0.1, r/rs = 0.1 — 1 and 
r/vs — 0.5 — 5, in order to probe different regimes of sampling 
density. We identify all the particles residing in these radial ranges 
and compute the ratio of the DTFE estimate of density to the NFW 
expected value for each particle in each radial bin. The results for 
the three different radial bins together with that for all the parti- 
cles within r-200 are shown in Figure |2] Despite the fact that the 
different radial bins have different sampling density, the one-point 
distributions in different radial bins are all the same and are consis- 
tent with that obtained for a Poisson point sampling of a uniform 
distribution. 

3.2.2 Dependence on mass resohition 

Halos identified in N-body simulations (or galaxy groups found in 
surveys) consist of different numbers of particles. In order to under- 
stand how the one-point distribution of DTFE reconstructed density 
fields depends on the number of particles used to resolve them, we 
generate a set of five NFW halos each with the same parameters 
M200 = 10«Mq, '"200 = 50kpc and c = 5, but with different 
numbers of particles: N20Q = lO'', 10^ lO'', 10^, and 10^, re- 
spectively. The one-point distribution function of the DTFE recon- 
structed density field of a N200 ~ 10^ halo will be noisier com- 
pared to a A'^200 ~ W^' halo due to effects of discreteness. In order 
to take the discreteness effects into account we generate different 
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Figure 4. The left and right panels show the density distributions constructed with DTFE for the MS and the MS-II, respectively. After computing the DTFE 
density for all the particles, we identify the particles residing in FoF halos of different mass ranges, and the particles which do not reside in any halo, and then 
compute the density distribution for each of the components separately. The distributions for non meinber particles and for particles in different halo mass bins 
are shown with different colours and symbols (as indicated in the panels) together with the distiibution for all the particles (black squares). 



numbers of NFW halos for each resolution: 1 with A^20o = 10 , 10 
with N200 = 10^ 10^ with iV200 = lO", 10^ with 7V200 = 10^ 
and 10^ with A^200 = lO'^ particles. One then has same total num- 
ber of density estimates for each resolution, allowing a straightfor- 
ward comparison of the one-point distribution functions at the same 
noise level. 



In the top panel of Figure [5] we show the one-point distribu- 
tion of p/ p-NFwij's) for all the particles within r2oo for NFW halos 
with A'^200 = 10^ and N2Q0 = 10^ particles. The distributions look 
very similar, except that the high density tail of the distribution for 
the N2m = 10 halo shows a slight shift towards lower density as 
compared to the N2(io = 10^ halo. 



We also compute the cumulative mass fractions as a func- 
tion of density from the distributions of the Ai'200 ~ 10^ and 
A*'2oo ~ 10® halos, and compare them in the bottom panel of Fig- 
ure [3] The theoretical prediction is shown with a solid line. This 
quantity gives the fraction of the virial mass contained when the 
NFW density profile is integrated from the center of the halo up to a 
certain density. The cumulative mass fraction as a function of radius 
is also shown, by indicating the radius along the x-axis on top. In- 
terestingly, the plot shows that even with 7V200 = 100 one can quite 
reliably reproduce the cumulative mass fraction from the measured 
one-point density distribution of the NFW halo. The = 10® 
halo overpredicts the cumulative mass at p > ps which roughly 
corresponds to r < rs/2. This overprediction is related to the large 
scatter in the DTFE density estimates which causes a broad range 
of radii in the halo to contribute to any individual density bin. Fur- 
thermore, the range of radii contributing to a density bin becomes 
even broader near the center of the halo due to the shallower profile 
of the inner region. The core of the iV2oo = 10® halos are resolved 
better relative to the halos sampled with A'^200 = 10^. The halos 
with A^200 ~ 10^ particles tend to underestimate the densities near 
the center. Here the boost due to scatter in densities is somewhat 
compensated by the poor resolution, enabling it to nicely but ar- 
guably misleadingly recover the analytic mass profile even near the 
center of the halo. 



Halo mass range Fraction (%) Fraction (%) 
(inh-iMe) in MS in MS-II 
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7.90 


1012 
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8.391 


1013 


- 10" 


5.981 


9.103 


10" 


-IQi^ 


7.639 


10.435 




-IQie 


9.972 





Table 2. Fraction of particles in each halo inass bin that account for the 
excess in the actually measured distribution function in the corresponding 
halo mass bin relative to the theoretical distribution functions. 

3.3 The one-point distribution function in the Millennium 
and Millennium-II simulations 

In the halo model, it is assumed that all the matter in the Universe 
resides in a halo of some mass. With this assumption, the whole 
matter distribution can be represented as a superposition of a set of 
halos in different mass ranges. To define the model one only needs 
to specify the density profiles of halos and the halo mass function. 
The density profile of dark matter halos can be described by the 
NFW profile (equation 0, which is a function of radius and mass 
of the halo. The concentration c = r2Qo/rs depends weakly on 
halo mass, and assuming the mass-concentration relation is known, 
one can write down the density p{r, M) for any particle residing at 
a radius r of a halo of mass M. For a smooth NFW halo the proba- 
bility distribution function P{p) is simply given by the fraction of 
the volume at density p, i.e. 

V dp r^go dp 

For an NFW halo, = ^^4^- So for a given mass of the halo 
and a specified mass concentration relation one can analytically cal- 
culate the density probability distribution function of the halo. 

We now contrast this expected density distribution with actual 
measurements for the MS and MS-II when the density field is con- 
structed with the DTFE. The probability distribution functions for 
all the dark matter particles at redshift z = from the MS and 
MS-II are shown as black curves in the two panels of Figure|4] The 
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Figure 5. Compaiison of DTFE density distributions for all the particles, all the non-halo particles, and all the particles in different halo mass bins in MS and 
MS II. The bottom four panels show different halo mass ranges, as labeled. The arrows in each panel show the density corresponding to the softening radius 
in MS and MS-II for each mass bin. The top left shows all the components together (i.e. all the particles in the simulations) while the top right compares only 
the non-halo particles. 



plots show that with the DTFE we are able to recover density values 
spanning about ~ 1 1 orders of magnitude, from underdense voids 
to the fully collapsed halos. The distribution is flat over nearly six 
orders of magnitude in density, and it does not exhibit any apparent 
linear to non-linear transition. The lower limits of the distribution 
functions shown in these figures have a Poisson error of ~ 3%. 

To isolate the contribution from collapsed halos we consider 
dark matter particles according to the masses of their host FoF ha- 
los. We identify the host halos by cross-matching particles IDs with 
FoF group IDs in the simulation. We also identify all the particles 
which are not part of any FoF halo with 20 or more particles. The 



density distributions for each of these cases are computed sepa- 
rately. The results for the different components are shown together 
in Figure |4] It should be noted that even though the low mass halos 
are more concentrated compared to their high mass counterparts, 
the distribution systematically extends to higher density values for 
the higher mass halos. This is opposite to what one would naively 
expect. This result could be explained by the fact that the smaller 
mass halos are sampled with fewer particles which poorly resolve 
the highly concentrated cores in these halos. In addition, the effect 
of gravitational softening (which introduces a soft core in the halo) 
in N-body simulations is expected to be more severe for less mas- 
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Figure 6. Comparison of our analytical halo model predictions for the one-point distribution of the density with the direct numerical simulation results summed 
over an equal range of halo mass bins (top row), or with the similarly summed Monte Carlo realizations of NFW halos over the same mass range (bottom 
row). For the MS we carry out the sum over 5 equally spaced different halo mass bins, each spanning a decade in the halo mass range lO" - lO^*^ h'^MQ. 
Similarly, for MS-II the sum is caiTied out over 6 equally spaced different halo mass bins in the halo mass range 10® — 10^^ h~^MQ. In both cases, we show 
results for the MS (left panels) and the MS-II (right panels) as red dashed lines, while the analytical halo model is shown with soHd black lines. 



sive halos as the softening length is a relatively larger fraction of 
their virial radii. The results for the MS-II in the right panel are 
quite similar, except for the fact that the overall shape of the distri- 
bution function at intermediate densities is somewhat different than 
that found for the MS. The distribution in this case is curvier and 
not quite as flat. This presumably reflects the larger fraction of par- 
ticles bound in halos in the MS-II compared with the MS. It is also 
noticeable that at a given halo mass the density distributions extend 
to high densities in MS-II. Again this generally reflects the larger 
softening of the MS. 

More detailed comparisons between MS and MS-II are shown 
in Figure|5] As the DTFE does not use any specific length scale for 
smoothing it is assured that the truncation in the tails of the distribu- 
tion is not a result of a spatial resolution limit imposed by the den- 
sity estimator. However, the intrinsically limited volume and mass 
resolution of the simulation is expec ted to introduce an undersam- 
pUng of the tails of the distribution ( IColombiiri994 iBagla & Ravi 
I2OO5I) . Smaller volumes undersample rare events in both overdense 
and underdense regions. A similar effect is caused by lower mass 
resolution primarily due to its limitation in resolving lower mass 
halos and smaller voids. MS-II has 125 times better mass resolu- 
tion and 5 times smaller softening length than MS enabling it to 
resolve the highly concentrated smaller mass subhalos. It can be 
seen in the top left panel of Figure |5] that the distributions in MS 
and MS-II are quite similar, apart from the fact that the distribu- 



tion for the MS-II extends slightly more on both the low density 
and high density ends. The effects of finite volume and finite mass 
resolutions somewhat compensate each other in the MS and MS-II 
but the slightly extended tail of the distribution at the both end in 
the MS-II suggests that the combined effects of finite mass resolu- 
tion and gravitational softening are more important than the effect 
of finite volume. 

The results for the non-halo particles are shown in the top right 
panel. The distribution in the MS-II has lower amplitude than the 
MS in intermediate density ranges because of its higher mass reso- 
lution and smaller softening, which allows it to resolve more bound 
objects at small mass. It is also interesting to note that the high den- 
sity tail of the distribution for the non-halo particles extends up to 
~ 10'' pc suggesting the presence of some high density sites even 
outside the bound FoF halos. The distributions in different individ- 
ual halo mass bins in MS and MS-II are shown in the two mid- 
dle and two bottom panels. The MS-II distributions have a sharper 
low-density cut-off and a more extended high density tail than in 
the MS in all cases. Any halo of a particular mass is sampled with 
125 times more particles in MS-II than in MS. This gives MS-II 
better power to resolve small subhalos and the boundary of the FoF 
groups. Further, the softening length in MS-II is 5 times smaller 
than in the MS reducing softening effects. These effects are more 
dominant in low mass halos and can be clearly seen as larger shifts 
in the lower halo mass bins. 
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Figure 7. The top left panel shows the one-point density distribution function (black squares) of MS particles residing in halos in the mass range 10^^ — 
10^^ h^^MQ. For comparison, the one-point distribution computed from Monte Carlo realizations of a NFW halo with mass 

10^2 /i-IMq is scaled 

according to the total number of particles present in halos in the mass range 10^^ — 10^^ h~^MQ. The distributions computed by restricting particles to 
within r G [0, rvir], [le, '"vir] and [3e, rvir] are shown in different colors and symbols. The other panels, except for the one in the bottom light, show the 
results for halo mass ranges 10^^ — 10^'^, 10^^ — 10^^, 10^'* — 10^^ and lO'^^ — 10^'' /i~^Mq, respectively. Finally, the bottom right panel shows the sum 
of scaled distributions of NFW halos in different mass bins, and adds it to the distribution of all the non-halo particles in the simulation. The full one-point 
distribution of all the particles in the MS is also plotted (black triangles), for compaiison. We note that while summing the distributions of Monte Carlo NFW 
halos over different halo mass bins we considered only particles within r G [le, Tvir] of the halo for the two lowest mass bins. For the rest of the bins all 
particles are considered. 



In the two top panels of Figure |6l we compare our analytical 
halo model for the one-point distributions of the density summed 
over different halo mass bins against that directly obtained for the 
MS and MS-II. For MS and MS-II, we sum the results for 5 and 6 
equally spaced different halo mass bins each spanning a decade in 
the halo mass ranges 10"-10^'^ /i~^Mq and 10^-10^^ h'^MQ, 



respectively. The excess in the observed distribution compared to 
the model prediction is higher in the MS-II than in the MS which is 
most likely related to the relative abundance of substructures. We 
explore this issue in detail in the remaining part of the paper One 
can also clearly see a larger suppression of the high density tail of 
the distribution in the MS compared to the MS-II due to its larger 
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Figure 8. The same as Figure|7] but for the MS-II. 



softening length. The shaip drop in the analytical predictions in the 
low density regime corresponds to a truncation of the halos at their 
virial radii. In the bottom two panels of Figure |6] we also compare 
the analytical predictions for the one-point distribution function of 
NFW halos summed across the different halo mass bins used in 
the analysis, and the combined one-point distribution function of 
DTFE densities computed from Monte Carlo NFW halos across the 
same halo mass bins. We use the same mass-concentration relation 
as employed for the analytic estimates. We see that the analytical 
one-point distribution function is quite well described by the results 
from Monte Carlo simulations, apart from the fact that the analyt- 
ical predictions show a sharp drop in the distribution function due 
to the truncation of all NFW halos at r2oo- The DTFE densities do 
not show this sharp drop due to the Poisson sampling of the ha- 



los. The slightly more extended high density tail seen in the MS-II 
comes from its ability to incorporate lower mass halos which are 
more concentrated. Thus, the DTFE traces the analytical one-point 
distribution function of the densities quite well, and the amount of 
excess (substructure) and the shape of the high density tail of the 
distribution in the N-body simulations are governed by the mass 
resolution, gravitational softening and simulated volume. 

We have also modeled the one-point distributions for different 
mass bins in the N-body simulations with ideal spherical NFW ha- 
los of similar masses. For this we generate 10'' , 10^, 10^ , 10 and 1 
Monte Carlo realizations of NFW halos with masses lO^^h-'^ Mq, 
lO"/i"^M0, lO"/i-^M0, 10^^ h-^ Mq and IO^^/j-^Mq, 
respectively. We use the mass-concentration relation given by 
iNeto et all bOOTh and the particle mass 8.61 x 10* Mq of the 
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MS. The concentration of halos depends very weakly on mass, and 
for simplicity we assume that the concentration does not change 
significantly within each of the narrow halo mass ranges. Different 
numbers of halos for the different mass ranges are constructed to 
account for the effects of discreteness i.e. to incorporate the fact that 
a 10"^ Mq halo is resolved with 10"^ ^"^ times more particles 
than a 10"^ Mq halo for any given values of ni and 712 (where 
ni > n2). We use all the density estimates from 10"^ ^"^ simulated 
NFW halos with mass 10"^ Mq to compute their one-point 
distribution, ensuring that the same number of density estimates are 
used to determine the one-point distribution of density for NFW 
halos with different masses. We followed the same approach for 
the MS-II as well, except that here we can go down another two 
decades in halo mass. We use the same mass-concentration relation 
for the MS-II, extrapolated to lower halo masses as needed. 

We compute in this way the one-point distribution of Monte 
Carlo NFW halos with different masses, and then scale them ac- 
cording to the total number of particles present in different halo 
mass bins as directly measured in the MS, in order to compare the 
theoretical predictions of this simple halo model with the measure- 
ments directly obtained from the MS (Figure |7). It should be noted 
that we have not used the halo mass function to model the one- 
point distribution here. Instead, we have used the total number of 
particles in a halo mass bin to predict the expected one-point distri- 
bution from halos in that bin. The results of this comparison for the 
MS and the MS-II are shown in Figures |7] and [8] respectively. 

We find that the one-point distribution in each mass bin is re- 
produced nicely at inteiTnediate densities when the results for the 
Monte Carlo NFW halos in that mass bin are scaled according to 
the total number of particles found in the N-body simulation in the 
same mass bin (Figures |7] and [§}. The low density part of the dis- 
tribution in each mass bin is however missed by the theoretical 
toy model. This is to be expected as we are only considering the 
particles within the virial radii of the spherical mock NFW halos, 
whereas in reality the FoF halos in N-body simulations extend be- 
yond their virial radii and generally have quite irregular shapes near 
their edges. At higher densities, the lower mass bins show higher 
values of the PDF for the theoretical one-point distribution com- 
pared with the N-body simulations. This is because the mock NFW 
halos do not involve any softening whereas a gravitational soften- 
ing is present in the N-body simulations from where the FoF ha- 
los are identified. In order to explicitly check for the impact of the 
softening length (e) in N-body simulations we have made an ex- 
periment where we prevented that particles in the mock halos are 
placed within le from their centers when computing the one-point 
distribution for all NFW halos in each mass bin. Interestingly, when 
we limit the particles to the radii outside of the softening range, the 
high density tail is nicely consistent between the N-body simula- 
tions and the mock halos. It thus seems clear that softening primar- 
ily affects the high density tail of the one-point distribution by sup- 
pressing the highest density values. This simple explanation does 
not work as well for the highest halo mass bins, presumably be- 
cause their structure is less strongly affected by the softening and 
they feature much more halo substructure. 

It is interesting to note that the one-point distribution in the 
higher halo mass bins shows an excess over the theoretical pre- 
diction. The fractions of particles in each halo mass bin which ac- 
count for this excess are different and increase with increasing halo 
masses. Specific numbers for our measurements are reported in Ta- 
ble 2. In a smooth NFW halo, there are no substructures whereas the 
halos foiTned in N-body simulations host numerous small subhalos 
that typically account for a few percent up to ~ 15% of the mass. 



The excess we find in the one-point distributions is most likely the 
direct consequence of the presence of substructures in the massive 
dark matter halos formed in N-body simulations. The excess ac- 
counts for up to ~ 8 — 10% of the total particles in the highest 
halo mass bin, consistent with typical substructure mass fractions. 
The excess decreases with halo mass and is completely absent or 
almost negligible in the lowest mass bins. Again, the decrease in 
the amount of substruct ures with decrea sing halo mass is consistent 
with previous findings jOao et ^l2004h and with expected numer- 
ical resolution limitations. Further, it is to be noted that the MS-II 
shows a higher substructure abundance compared to the MS in each 
halo mass bin. This is related both to the higher mass resolution and 
to the smaller softening length in the MS II which enables it to re- 
solve lower mass subhalos. 

The effect of substructures on the one-point distribution is 
shown explicitly in Figure |9] which focuses on a well resolved 
halo from the MS-II. The substructures i n the MS-II halo are here 
identified with the SUBFIND algorithm jSpringel et"ani200lh and 
are available from the MS-II database. We identified all the par- 
ticles within 1.5 r^ir around the center of that halo and computed 
the one-point distribution of the density. We then removed all the 
substructure particles within that radius and computed the distribu- 
tion function again. The results are compared in Figure |9] which 
clearly highlights the excess due to the presence of substructures in 
the halo. The substructures constitute ~ 6.5% of the total particles 
in this halo within our chosen radius, and this fraction is consistent 
with the values listed in Table 2. 

One should also keep in mind that the dark matter halos in 
simulations are not exactly spherical. Rather their shapes resemble 
in ge neral triaxial ellipsoids with a preference for prolate configura- 
tions ijDubinski & Carlberg |l99il: ICole & Lacevll996t|jing & Sutd 
l2002l : Vogelsberger et al. 2009). The shapes of the majority of 
the halos forming in N-body simulations are c haracterized by a 
mean axis ratio of about 1 : 0.74 : 0.64 |Kasu n & Evrard 200^; 
iBailin & Steinmetzl2005l : lAllgood et al.l2006h . Fitting a NFW pro- 
file to such halos inevitably involves spherical averaging. The 
spherical averaging of a triaxial halo will introduce systematic dif- 
ferences in the density estimates in different parts of the halo. Some 
parts of the halo have in reality densities which are larger/ smaller 
than the spherically averaged densities we try to reproduce in our 
mock models. This deviation in the densities could be quite high 
depending on the triaxiality of the dark matter halos and it is im- 
portant to investigate if the excess in the observed one-point dis- 
tribution of density is related to this issue. In order to test this we 
simulated a triaxial NFW halo of mass M200 = 10^'' h-^ Mq with 
axis ratios a : b : c = 1.5 : 0.888 : 0.75. To mimic the effect of 
sphericalization we randomize the azimuthal co-ordinates (6, (j>) of 
all the particles in the halo. The one-point distribution function of 
the triaxial halo before and after sphericalization are compared in 
Figure [Tolwhere we can see that this effect does not make a signifi- 
cant change in the one-point distribution function of density despite 
our choice of a deliberately extreme (yet still possible) axis ratio. 
This result suggests that the observed excesses in the one-point dis- 
tributions (Figures |7] and [§} are a direct outcome of the substruc- 
tures present in the dark matter halos, and that the halo shape plays 
only a very subdominant role. 

As an aside, we note that the dependence of substructure abun- 
dance on mass resolution has important implications for predic- 
tions of the expected extragalactic gamma ray background hypo- 
thetically caused by the self-annihilation of dark matter particles in 
dark matter halos. The rate of WIMP annihilation and hence the 
intensity of the annihilation radiation in a dark matter halo is pro- 
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Figure 9. The one point distribution of density for a lialo from tlie MS-II for 
all the particles residing within 1.5r200 around the halo centre. For com- 
parison, the red dashed line shows the one-point distribution of the same 
halo after all its substructures within 1 . 5 r200 has been removed. The virial 
radii of this halo is 0.60 Mpc, and a total of 8079161 particles reside 
within 1.5 r'200. out of which 528753 (~ 6.5%) are subhalo members. 




1 10 100 1000 10"* 10^ 10^ lo^ lo" io^ 



P/Pc 

Figure 10. One-point density distribution function of a Monte Carlo real- 
ization of a triaxial NFW halo (black) of mass A/200 = 10^'* hT^Mqy, 
represented with particle mass 8.61 X 10^ h~^MQ. The density distii- 
bution function of the same halo after spheiicalization is shown in red. To 
avoid any boundary effects only the particles inside r2oo ^re used in both 
cases. 



portional to the volume integral of the square of the dark-matter 
density. The existence of lumps or subhalos in dark matter halos is 
expected to enhance the ann ihilation rate by a signifi cant boost fac- 
tor ^ Bergstrom et al. 199 ^ IStoehr et al.) 12003: Kou shiappas et al.l 
2004; Fieri etal. 2005; Kuhlen et al."2008; Spring el et alj l2008l : 
Kuhlen et al. 2009; Kamionkowski et al. 201C). The smallest sub- 
halos are the densest sites in the halo and the predicted boost factor 
depends significantly on how well these subhalos are resolved in 
the N-body simulations. Given that halos in the MS-II host more 
resolved substructures than the MS indicates that the theoretical 
prediction of the intensity of annihilation signal from any simu- 
lated dark matter halo depends on the resolving power imposed 
by the softening length and by the finite mass resolution. One has 
to correct for this by extrapolating the relations between substruc- 
ture abundance and mass resolution (at a fixed softening length) to 
lower particle masses in order to make a prediction about the signal 
expected in reality. 

We combine the results of our Monte Carlo simulations of 
NFW halos in different mass bins, and add also the one-point dis- 
tribution of the non-halo particles found in simulations to obtain a 
halo model for the full one-point distributions in the MS and MS- 
II. The results are shown in the bottom right panel of Figures |7]and 
[S] The dip in the middle of the distributions results from the trun- 
cated spherical boundaries of the idealized NFW halos, which does 
not take into account the fact that the real halos in simulations are 
far more irregular at their boundaries and extend to lower densi- 
ties as well. At higher densities, there are nevertheless also differ- 
ences between the direct results of the N-body simulations and the 
Monte Carlo model. Clearly, larger differences are seen in the case 
of the MS-II compared with the MS, an effect that we attribute to 
the higher abundance of substructures in the MS-II. 

Finally, we compare the distribution of DTFE and SPH densi- 
ties from the MS and MS-II in Figure[TT]to check how the choice 
of density estimator influences the tails of the distribution. The 
smoothing lengths in SPH are chosen such that the sum involves 
32 nearest neighbours. This figure shows that for the lowest density 
regions or voids the DTFE and SPH smoothing give very similar re- 
sults. At very low densities, the slightly higher values found in the 
DTFE distribution come probably from a more accurate representa- 



tion of void boundaries than in SPH. Similarly, the extended kernel 
of SPH smoothes out high density regions like filaments or halo 
centers leading to an underestimate of the PDF in such regions. In 
both the MS and MS-II there is a small bump in the density range 
0.1 — 5pc displaying an excess of the SPH density estimate as com- 
pared to the DTFE density measurement. 

In Figure [12] we show how the DTFE and SPH density esti- 
mates compare with each other at different densities in the MS- 
II. It can be clearly seen that for most of the particles, DTFE 
and SPH give comparable results, while residual systematic dif- 
ferences show up as oppositely skewed tails at high and low densi- 
ties. These systematic differences in low and high density regions 
account for the differences in the SPH and DTFE density distri- 
butions. A detail comparison between DTFE and SPH is given in 
|Pelupessv]e t al. (2003 ) and our findings are consistent with theirs. 
Recently, Abel et al. (2012) pointed out that even the VTFE un- 
derestimates the densities in regions around filaments and sheets, 
as compared to a novel technique that more accurately represents 
dark phase-space sheets. So besides factors like mass resolution, 
gravitational softening and simulated volume, the choice of den- 
sity estimator plays a crucial role in shaping the tails of the density 
distribution. 



4 CONCLUSIONS 

The present analysis shows that the part of the one-point distribu- 
tion function represented by collapsed halos can be quite well de- 
scribed by a simple superposition of a set of NFW halos over dif- 
ferent mass ranges. However, the one-point distribution functions 
in N-body simulations also show a prominent hump when individ- 
ual halo mass bins are considered, especially for the more mas- 
sive halos. This excess with respect to the distribution obtained for 
smooth NFW halos originates in the substructures present in the 
massive dark matter halos. The amount of resolved substructures 
depends on the mass of the halo, and especially on the finite mass 
resolution of the N-body simulation. Further, the gravitational soft- 
ening suppresses the high density tail of the one-point distribution 
in halos, introducing a soft core which is more noticeable in smaller 
mass halos. Both of these effects imply that the high-density tail is 
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Figure 11. Comparison of the density distributions measured with DTFE 
and SPH smoothing. The top and bottom panels show results for the MS 
and the MS-II, respectively. 



still underestimated both in the direct N-body simulations and the 
analytical halo model. 

We find that finite simulation volume, finite mass resolution, 
gravitational softening as well as the method for estimating the den- 
sity field all influence the tails of the measured one-point density 
distribution. We note that this distribution function is a particularly 
important simulation prediction, as it, for example, determines the 
intensity of the WIMP annihilation signal from a representative vol- 
ume, which sensitively depends on the ability to resolve the abun- 
dant yet dense small-mass structures. Our analysis with the DTFE 
in the MS and MS-II suggests that the effect of finite mass resolu- 
tion and gravitational softening are the primary limitations rather 
than a finite simulation volume. Also, it is worthwhile to employ 
the DTFE techniques instead of simpler schemes for density recon- 
struction such as SPH-like smoothing, due to its sharper resolving 
power. 
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Figure 12. Two-dimensional histogram of the ratio of DTFE to SPH den- 
sity estimated for the MS-II, as a function of the DTFE density. A total of 
100 X 100 pixels are used to estimate the histogram. Different contours 
correspond to the pixels with a count equal to the numbers used for la- 
beling the respective contours. Different coloured regions mark the pixels 
with counts intermediate to their bounding contours. The numbers shown 
inside the coloured regions give the total cumulative counts in all the pixels 
contained inside the outer contours of the corresponding regions. This plot 
clearly shows the systematic overestimate and underestimate of densities by 
SPH over DTFE at low and high densities, respectively. 
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